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Abstract: We study the stationary nonequilibrium states of N point particles 
moving under the influence of an electric field E among fixed obstacles (discs) 
in a two dimensional torus. The total kinetic energy of the system is kept 
constant through a Gaussian thermostat which produces a velocity dependent 
mean field interaction between the particles. The current and the particle dis- 
tribution functions are obtained numerically and compared for small |E| with 
analytic solutions of a Boltzmann type equation obtained by treating the colli- 
sions with the obstacles as random independent scatterings. The agreement is 
surprisingly good for both small and large N . The latter system in turn agrees 
with a self consistent one particle evolution expected to hold in the N — > oo 
limit. 
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1. Introduction 

In this note we continue our study of the stationary nonequilibrium states (SNS) of 
current carrying thermostatted systems. In part I [1] we described extensive numerical 
and analytical investigations of the dependence of the current on the electric field for a 
model single particle system introduced in [2] and previously studied in [3] . Here we study 
a generalization of that model to N particles introduced in [4]. The particles, which have 
unit mass, move among a fixed periodic array of discs in a two dimensional square A with 
periodic boundary conditions, see Fig. 1. They are acted on by an external (electric) field 
E parallel to the x-axis and by a "Gaussian thermostat". (The discs are located so that 
there is a finite horizon, i.e. there is a maximum distance a particle can move before hitting 
a disc or obstacle). 



Fig. 1: General billiard structure with discs of radius Ri and R2 in a periodic box 
with side length 2L, N = 3 particles are shown. 



E 




2L 



2 



The equations of motion describing the time evolution of the positions and velocities 
Vj, i = 1, N, are: 

cu =Vi qi = (qi, x , q%, y ) £ A' 

(1.1) 

v, =E - a(J, Lf)v; + F o6s (q,) 
where 

T F 1 N 1 N 

i=l i=l 

Here A' = A\P, with £> the region occupied by the discs (obstacles) and F obs represents 
the elastic scattering which takes place at the surface of the obstacles. The purpose of 
the Gaussian thermostat, represented by the term a (J, U)v in eq.(l.l), is to maintain the 
total kinetic energy 1/2 J2iLi v f constant, i.e. U = Vq. It also has the effect of making the 
flow $t generated by eq.(l.l) on the (4iV — 1) dimensional energy surface non Hamiltonian 
when E^O. In fact the phase space volume contraction rate is given by cr(X) = — (2iV — 
l)oi(J, U). Another effect of the thermostat is to effectively couple all the particles in a 
mean field way, a(J,U), depending only on the total momentum of the particles. Note 
that this is the only coupling between the particles in this system. 

The change of variables, q^ — > q^/L, Vj — > Vj/^o, t — > tv /L and E — > EL/i^, where 
2L is the length of the box, leaves eq.(l.l) unchanged, so that the motion of the system 
takes place on Sn = (A') N x S N , where S N = {v,| J^Zi v f = N }- We shall denote by 
X e Sn a point in the phase space of the system. In these units we took U = 1, R\ = 0.39, 
R 2 = 0.79, and A is thetorus of side 2 1 . 

Our main interest is in the SNS of this model system. To be more precise let no(dX, N) = 
Po(X; N)dX be an initial measure symmetric in the {qi,Vj} and absolutely continuous 
1 See [1] for an explanation of these values. 
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with respect to the Liouville volume dX projected on Sn- The time evolved measure 
Ht{dX, E; N) is still absolutely continuous with respect to the Liouville measure with den- 
sity pt(X, E; AT) for any fixed time t. The SNS is expected to be described by an SRB 
measure fi + (dX, E; N) , given by the weak limit, as t — ► oo, of /j, t (dX, E; N) , when it 
exists. This limit measure is in general not absolutely continuous with respect to the Li- 
ouville measure, due to the phase space volume contraction [5] , [6] . The existence of such 
a limit was proven, for A" = 1 and |E| e [0, Eq] (Eq small) in [3], but no such result is 
available for A" > 2, because of the lack of uniform hyperbolicity for the zero field system. 
On the other hand our computer simulations of the dynamics, for A" ranging from 1 to 
50 and E from 0.04 to 1.0, strongly support the belief that there exists a unique limiting 
measure [i + (dX, E; N) up to quite large values of |E|, say |E| = E < 1. We expect however 
that the projection of /i + (dX, E; N) on the one particle phase space A' x fi(jv), where fi(jv) 
is the ball |v| < VN, will yield a one particle density / + (q, v, E; N) absolutely continuous 
with respect to dqdv; this is proven, for instance, for coupled Arnold's cat maps [7]). 

To obtain information about / + we considered first the case of weak fields. It is tempting 
to think that for E — > the singular set on which fx + is concentrated will be spread out 
more or less uniformly on Sn so that fi + will approach weakly the microcanonical measure 
on the energy surface Sn- this measure is certainly invariant for the dynamics at E = 0. 
If this were the case then / + (q, v, E; N) would approach, as E — > 0, the equilibrium one 
particle density obtained from the projection of the microcanonical measure: for large A" 
this would be close to the Maxwellian distribution with unit variance 2 . We ran computer 
simulations for values of the field between 0.04 and 0.12 and N = 2,5 and 50. In all cases 

2 Note that for large N the Maxwell distribution is typical for points on the energy 
surface, i.e. the set B on Sn for which /+ is not a Maxwellian has measure (w.r.t. 
dX). Of course since /i+ is singular w.r.t. dX this need not to be the case here. 
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we found a one particle distribution that is far from the projection of the micro canonical 
distribution. Furthermore this distribution appeared to have only very slight dependence 
on E for those values of the field; so it appears that there is a well defined limit of 
/ + (q, v, E; N) as E — > 0, and that this limit is not the projection of the microcanonical 
measure: there are correlations between the velocities of the particles induced by the field, 
beyond those corresponding to the energy constraint, which remain when E — > 0. 

This deviation from the microcanonical distribution is reflected also in the behavior of the 
average current per particle in the steady state, given by j(E, N) = J v/ + (q, v, E; N)dqdv 
as E — > 0. We studied j(E, N) numerically as a function of E and N, see Fig. 2 and Fig. 3. 
In the following we will always assume that the electric field is along the positive x-axis, 
E = El x . This implies that the y component of j(E, N) is zero for symmetry reason. We 
will denote the x component of the current by j(E, N) and call k(E, N) = j(E, N)/E the 
conductivity. The dependence on N for E — > should be given by the Green- Kubo formula 
for the zero field conductivity when the dynamics of the particles are independent. A 
straightforward computation then shows that the zero field conductivity of the iV particles 
is: 



k(0,N) = CW(0)k(0,1) 



(1.3) 



with k(0, 1) given by the diffusion constant of Bunimovich and Sinai [8] and 




(1.4) 



For the microcanonical distribution we easily find: 




(1.5) 
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which is inconsistent with our data although the form of the dependence on N appear to 
be similar, see sect. 2.1. 

Let us consider now the behavior of our model system in the limit N — > oo. As the 
particles interact only through their average velocity 3(X(t)) it seems reasonable to expect 
that, for N — > oo, J will stop fluctuating, i.e. that for "well behaved" initial distributions 
[9], [10], [11] 

J(X(f))— >j t = J v/ t (v,E)dv (1.6) 

where / t (v, E) = lirn/v >oc /t(v, E; AT). If this were true in a sufficiently strong sense it 

would lead to an autonomous Vlasov type equation [9], [10], [11] for f t where v would be 
computed self consistently from the (irreversible) dynamics 3 



v = E-A(t)v + F 0&s (q) (1.7) 

with X(t) = E • } t . The difficulty with proving this behavior, as compared to the [9] case, 
is that trajectory X(t) and thus also 3(X(t)) is not smooth for finite t. The problems are 
compounded when we consider the t — ► oo limit corresponding to the SNS. 
Based on numerical evidence we nevertheless believe that 

lim /+(v,E; iV) = /+(v,E) = lim / t (v,E) (1.8) 

where /t(v, E) is the solution of the Vlasov equation with a force given by the right hand 
side of (1.7), and we define for a given function g 



#( v ) = / #(q> v ) rf q • 

J A' 



3 The dynamics (1.1) is reversible in the sense that if TtX is a solution then TfRTfX = 
RX, where R reverses all velocities. 
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The integration over q is necessary, or at least desirable, since we expect the t — > oo limit 
of /t(q, v, E) to be singular with respect to dqdv as is the N = 1 reversible system (1.1). 
Its projection on the velocity is however expected to be absolutely continuous with respect 
to dv [3] [7] . Eq(1.8) is thus a form of the law of large numbers which should hold for 
smooth p (X,~E; N). Something like this was in fact proven by Ruelle for the stationary 
state under some hypotheses on the thermostatted dynamics [12]. To make contact with 
Ruelle's theorem it is convenient to think of A at as a torus of length 2LN along the y-axis 
(perpendicular to E) and length 2L along the x-axis. This does not change the dynamics. 

To get some analytical handle on the form of the reduced distributions in the SNS we 
investigated a model system in which the deterministic collisions with the obstacles are 
replaced by a stochastic process in which particle velocities get their orientations changed at 
random times, independent for each particle. This yields a Markov process which replaces 
the continuity equation for p t (X,E;N) by a linear Boltzmann-like equation, see [13]. We 
can write either of these equations in the symbolic form: 



|o<(Q, V) + £ A {v,p t (Q, V)} + £ °- { [E - a(V)v,] p t (Q, V)} = (|^ , 

(1.9) 

where we have set X = (Q, V) (and dropped the explicit dependence on E and AT). The 
term on the right hand ( % L ) represents either the effect of deterministic collisions with 

V / coll 

the obstacles as given by (1.1) or a collision operator independent of Q, see (3.1). A similar 
ansatz for the irreversible dynamics (1.7) leads to a Boltzmann-Vlasov equation for the 
one particle distribution. These equations can be solved analytically as a power series in 
E and/or numerically. This is described in sect. 3. 

In sect. 4 we compare some of the moments, including the current, of the determinis- 
tic distribution / + (v, E; AT) with those of the stochastic one. We find surprisingly good 
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agreement once the mean free path appearing in the Boltzmann-like equations is properly 
interpreted, see sec. 4.2. We note however that a direct computation of the distribution of 
free paths in the dynamical system (1.1) shows that it is far from being exponential, which 
is the basic assumption of the Markov process. We therefore have no real explanation for 
the observed good agreement. We only note that some features of the stationary state 
appear rather robust with respect to the collision processes with the "obstacles" , yielding 
similar results for different distributions for the free path. In sect. 5 we discuss some gen- 
eral questions about the relation between this thermostatted model and the Drude model 
of electrical conduction in metals [14] . 

2. Numerical results 

Eq. (1.1) can be solved in terms of quadratures between collisions with the obstacles so the 
simulation consists mainly in computing the times of successive collisions. At each collision 
there is an instantaneous change in the velocity of the colliding particle and consequently 
also in the current J and thus in the thermostatted force acting on each particle. Assuming 
that the system is ergodic we can obtain information about the SNS from time averages 
over a single trajectory. In practice we used a few initial states and found a behavior 
consistent with this assumption. The relative simplicity of the dynamics enabled us to get 
fairly accurate results even for 50 particles with relatively small computing power. Our 
simulation were carried out on a Pentium PC. Error bars are computed by doubling the 
range of the fluctuations of the time average over the interval [0.9T, T] where T is the 
total number of collisions computed. After the change of variables described after (1.2) all 
quantities appearing in the graphs are adimensional. 

2.1. The current 

Let j(E, N) be the average current in the steady state 
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j(E, TV) = (J) M+ = J v/+(v,E, iV)dv . 



(2.1) 



with J defined in (1.2). As already noted, in all our computations the electric field is along 
the positive x-axis, E = El x , all densities are normalized and j(E, N) is the x-component 
of the current defined in eq.(2.1). 
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Fig. 2: Conductivity k(E, N) as a function of E for different N. 



In Fig. 2 we plot the conductivity n(E,N) = j(E,N)/E as a function of the field for 
different numbers of particles, A/=l,2,10,15,20,30 and 50. The averages were computed by 
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Fig. 3: k(E,N) as a function of iV _1 for different £7. Also plotted is the conductiv- 
ity obtained from eq.(1.3) using the actual distribution function, see next section, for 
E = 0.04 and compared with the value obtained by a direct simulation at the same 
field. Finally the highest line represents the conductivity obtained from eq.(1.3) using a 
microcanonical hypothesis. 



running simulations in which the total number of collisions with the obstacles varied from 
10 9 for AT = 1 to 10 8 for N = 50. 

We note that for very small fields the interaction among the particles is very small so 
that the invariant distribution is reached only after a very long transient time. 

Furthermore, although the current goes to as E — > 0, the fluctuations in the current 
are almost independent of E so that longer and longer simulations are required in order 
to distinguish the average from the fluctuations when E — > 0. For N = 2, 5 and 10 we 
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checked whether — - — > as E — > 0, as required by the symmetry of the problem if 
k{E, A) is differentiable at 0. While the results are not definitive they are consistent with 
such behavior. 

In Fig. 3 we plot the conductivity as a function of 1/A for a few selected values of 
the field. As can be seen there the behavior of k(E, A) can be well fitted for A > 2 by 
the following formula which is the analogous of eq. (1.3) with CV(0) given by (1.5) for 
E 7^ 0: k{E,N) = k(E) + c/N with k(E) = limAr^oo k(E,N) and c independent from 
E, at least within the accuracy of our computation. (The value of k(E, 1) is about 15- 
20% lower than that given by the formula, depending on E). For E = 0.04 we have the 
value of the conductivity for A = 2, 5 and 50 as well as the distribution / + (v, E; A). We 
can therefore check directly eq.(1.4) for E ^ 0. Fig. 3 contains both the values obtained 
directly and those obtained from eq.(1.4) for E = 0.04. The agreement is clearly very good. 
Finally plotted in Fig. 3 is the value of the conductivity at zero field obtained from eq.(1.5), 
i.e. assuming that the invariant distribution is microcanonical. Although this assumption 
is inconsistent with the actual numerical data, the behavior is qualitatively similar. 

The smoothness, or rather the lack of smoothness, of the current as a function of E 
for A = 1 was extensively discussed in [1] and related there to the discontinuities of the 
collision map. The data we have for A > 2 are insufficient to address this question. 
However it is expected that the stationary current will be smoother than it is in the one 
particle case, since it is averaged over all particles. 

2.2. Distribution functions 

To study the space independent part of the one particle density function, / + (v, E; A), 
it is convenient to switch to the variables r = |v| e [0,y/~N] and 8 e [— n, n] the angle 
between the velocity v and the x-axis. Expanding / + (v, E, A) in a Fourier series in 9, we 
have 
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f+(v,E;N) = Y,Mr,E; N) cos k9 , 



(2.2) 



k=0 



where only terms in cos kO appear due to the symmetry of the problem. Note that 
2TTripo(r, E; N) is the stationary probability density for the modulus of v while 



j(E,N) = 7T dr^^r^E-N) 
Jo 



(2.3) 
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Fig. 4: Plot of 2nripo (r, 2) for different values of £7. The straight dashed line is 
obtained from the microcanonical distribution, Eq.(2.4). The dotted line gives the result 
for the stochastic model 
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Fig. 5: Plot of 7rripi(r,E;2)/E for different values of E. The dotted line gives the 
result for the stochastic model 

In Fig. 4 we plot 2irrip (r, E; 2) for E = 0.04,0.08,0.12 while Fig. 5 is a plot of 
7r7"i/>i(r, E;2)/E for the same values of the field. Both appear to be almost independent of 
E for those values of E so we believe that Figs. 4 and 5 represent a good approximation 
for the limiting behavior E — > 0. Observe that, due to the symmetry E — > — E we expect 
the corrections to these functions to be of 0(E 2 ). For comparison we also plotted there 
the results obtained analytically from the stochastic model discussed in the Introduction 
and in section 3. 

In Fig. 4 we also plot the " micro canonical" density of |vi| obtained from the micro- 
canonical ensemble of 2 particles with vf +v% = 2. The microcanonical one particle density 
/micro( v ) is of course isotropic and the speed distribution, 27r|vi|/ m (|vi|, E = 0;2), is 
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27r|vi|/ m (|vi|, J E = 0;2) = i| Vl | J 5{y\ + v 2 2 - 2)rfv 2 = | Vl | W(2 - v?) , 



(2.4) 



where TC(x) is the Heaviside function. This is seen to be very different from what we 
obtain from our simulations or analytically from the stochastic model for E — > 0. We did 
a similar analysis for N > 2 and in Figs. 6 and 7 we present the corresponding results for 
N = 50. 
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Fig. 6: Plot of 2irripo(r, E; 50) for i? = 0.04. Also shown are the results from sim- 
ulations of (1.7) and from analytic solutions of the corresponding stochastic equation, 
Eq.(3.10). For comparison we also show the microcanonical result, corresponding to a 
Maxwellian. 
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Fig. 7: Plot of 7vri/ji(r, E; 50) /E and comparison with stochastic irreversible dynamics 
for E = 0.08 

2.3. The N = oo /zmzt 

As discussed in sec. 2.1, k(E,N) — > k(-E) as — > oo. We compared the k(E') obtained 
from our simulation, see Fig. 3, with that obtained from the irreversible eq.(1.7). A way 
to do this self- consistently would be to choose the parameter A in eq.(1.7) such that 



U(E) = J , 



dv\v\ 2 f + (v,E) = 1 

and show that for this value of A the conductivity k(E) for the system described by eq.(1.7) 
is equal to k(E). Rather than doing this, we took the k(E) deduced from the simulations 
as in Fig. 3 and used it to determine A, i.e. we set A = k(E)E 2 in eq.(1.7). We then 
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computed, via simulation of eq.(1.7), a new conductivity k(E). In Fig. 8 we compare k(E) 
and k(E). The agreement is very good. We observe that it follows from eq.(1.7)that 
E 2 k(E)/U(E) = A so that this agreement also confirms the self-consistency discussed 
above. 



y 



0.285 



0.28 



0.275 



0.27 



0.265 



0.26 



0.255 



0.25 



irreversible dynamics >- 
large N limit of reversible dynamics t - 



0.2 



0.3 



0.4 



0.5 



0.6 
E 



0.7 



0.8 



0.9 



Fig. 8: Comparison between the limiting value of the conductivity k(E) in the reversible 
model and in the irreversible model Koo(E). 



As for the reversible dynamics we can write 

oo 

f+(v,E) = J2 Mr, E) cos k6 (2.5) 

k=0 

In Figs. 6 and 7 we compare Itxt^q (r, E; 50) and irripi(r, E; 50) with 2irr<f>o(r, E) and 
nr(f)i(r, E) respectively. The agreement is very good. As we did for = 2 in Fig. 4 
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and Fig. 5 we also plotted in Figs. 6 and 7 the results obtained analytically from the 
stochastic model discussed in the Introduction and in section 3. In Fig. 6 we also plot the 
micro canonical density, i.e. a Maxwellian with < vf >= 1. 

3. Thermostatted Stochastic Evolution 

We now describe more precisely the stochastic model system in which the collisions 
between particles and obstacles are replaced by independent random scattering events. 
The model is specified by writing the right hand side of eq. (1.9), the evolution equation 
for the iV-particle phase space density of our system, which we now call F t (Q, V), to 
distinguish it from the mechanical pt(Q, V), as 

W^T) =rl £/ K_£) (F(Q , V ;,E)-F(Q,V,E)) d „ (3.1) 

V ot y coll ~-tAn-v,)<o 2 

In (3.1) n is a unit vector in the direction of the momentum transfer in a "collision", 
|n| = 1, v' = v — 2n(n • v^) and is identical to except for its i-th component which 
is replaced by v^. The coefficient / _1 multiplying the collision term is the inverse of the 
mean free path between collisions, a parameter to be specified. 

Eq. (1.9) together with (3.1) describes a Markov process in which particles change the 
directions of their velocities as if they were undergoing independent random collisions 
with "phantom obstacles" at a rate equal to / _1 |v| with a uniformly distributed impact 
parameter [15]. Between collisions the particles move according to eq.(l.l). This model 
can be thought of as, and presumably even proven to be, the Boltzmann-Grad limit of our 
system: i.e. , we place discs of radius R randomly in a square of side L with density p and 
then take R — > 0, p — > oo such that / = stays constant, see [16]. 

This system will, like our mechanical system, eq.(l.l), conserve energy, so setting J2 v i = 
N the evolution takes place on Sn- By general arguments [17], [18] we expect that this 
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system will, for F 7^ approach, as t — > 00, a unique stationary density F(V, E; AT) which 
will satisfy the equation 



f;^ {[E _ E . Jv , ]F( v, E;J v) }= (^™) 

i=l 1 ^ 



) 



coll 



(3.2) 



For small E we expand F(V, E; N) as a formal power series in E: 



CO 



F(V, E; JV) = F(R, 6) = # n ^ (n) (R, ©) 



(3.3) 



n=0 



where we have set Vj = (r; cos6>;, sin6>;) and R = (ri, . . . , rjv), Zli r i 2 = AT, © = 
(Oi, . . . ,$n)- Observe that in this way we get a singular perturbation problem because F 
multiplies the highest order derivative in eq.(3.2). Moreover F + (V, E; N) clearly depends 
only on E/l so that we can, for the time being, set / = 1. Finally we can write, as in the 
previous section, 



where we have again used the symmetries of the problem. 

Substituting (3.4) into (3.3) one gets a hierarchy of equations linking F^ (R, k) to 
F( n_1 )(R, k l ) where k l = (ki,...,ki + l,...,fcjv). From this, and from the fact that 
the kernel of the collision operator depends only on R we get that F^ n ^(R, k) = if 
|k| > n. F(°)(R, 0) satisfies the relation: 



N 




(3.4) 



d_ 

dr 



i 




(3.5) 



while for F^^R, l ) we get the equation 



^{{-v-^ FW{n ' ' ) + i FW{K ' oi) } 



= 



(3.6) 
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with U = Ylirf- Equations (3.5) and (3.6) are easily solved and, together with the fact 
that FW(R, 0) = give us F(R, 0) to first order in E 



N 



F(R,0) = C8(£r*-N) 



i=l 



3(2N-1)E ncosOj 2 

2N-1 ' a 2N + 2 ' ^ V ) 



(3.7) 



where C is a normalization constant. It is possible to write out the full hierarchy of 
equations for F^ n )(R, k) and see that they can be solved iteratively but it is not clear that 
this is useful. We shall therefore use eq.(3.7) to compare with our numerical data for small 
values of E. To do so we define the one particle distribution /(v, E; N) and develop it in 
a Fourier series exactly as in eq.(2.2): 

/oo 
dv 2 • • • dv N F(V, E;N) = J2 Mr, E; N)cos(ki0i) (3-8) 

k=0 

Before doing any comparisons we consider the stochastic version of the fi(v, E) obtained 
from the irreversible dynamics defined by eq.(1.7). Putting A = E 2 v, v to be set to R(E) 
when compared with the deterministic model, we get 

A{[E-,Vv]/.(v,E,Wte) p9) 

\ / coll 

where the collision term is again given by eq.(3.1) with N = 1. Observe that although 
eq.(3.9) contains three parameter {E, v and /) it depends only on El and vl~ x . Developing 
/i(v, E) in a power series in E we obtain in analogy to (3.7) 

/i(v,E) = Ce~^ ur3 (1 + 2uEr cos 6) + 0(E 2 ) (3.10) 

where C is a normalization constant. 

To compare fi(v, E) with the large iV limit of /(v, E; N) given in (3.8) and (3.9) we need 
to fix the parameter v (setting 1 = 1). This can be done self-consistently requiring that: 
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j |v| 2 /,(v,E)dv = l (3.11) 
Solving eq.(3.11) for v and using it to compute fi we expect that: 

lim f(v,E;N) = f i (v,E) (3.12) 

While we have not proven this equivalence we believe that it should follow from general 
considerations: it would follow formally from showing that, in the limit N — > oo, F(v, E; N) 
factorizes, as is usually the case for systems with mean field type interactions. This is 
certainly consistent with our numerical results. 

4. Comparison between the deterministic and stochastic time evolution 



4-1. The distribution of the modulus of v 

For N = 1 the exact solution, for E = 0, of both the stochastic and mechanical models 
is /(v, 0; 1) = 5(v 2 — 1). For N = 2, we are able to compute the one particle distribution 
from eq.(3.7). This yields 

Cr 

^o(r, E- 2) = r3 + (2 _ r2)3/2 + Om (4.1) 

where C is a normalization constant. This is plotted in Fig. 4 and one can easily see that 
the agreement with the numerical solution of the deterministic model is very good. 

A similar agreement is obtained for iV = 5 although, as already said we were not able to 
integrate eq.(3.8) for iV > 2 so that we computed this integral numerically by simulating 
the process associated to eq.(1.9) with collision term given by eq.(3.1). 
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Finally for N = 50 we see in Fig. 6 that our deterministic (1.1), stochastic (3.9) and 
irreversible (1.7) models give indistinguishable results. This certainly suggests the validity 
of (1.7) and (3.12) for large N. 



4-2. The first Fourier component of the distribution of v 

The analysis of the first Fourier component of the distribution of v is less straightforward 
because we must fit the parameter / appearing in eq.(3.1). In the stochastic system / 
represent the mean free flight of a particle. The concept of mean free flight is not uniquely 
defined for the mechanical model. For this reason we used / as a fitting parameter for 
matching V>i(r, E; N) with ipi(r, E; N). We will go back to the mechanical meaning of this 
parameter in the following section. The case N = 2 is reported in Fig. 5 where, for the 
periodic case, we used a field E = 0.04 and for the stochastic one we have the expression 

1 OR/ /7 r 2 

*<^ ;2 »-V(^bi^ +0(£) (42) 

with C the same costant appearing in eq.(4.1) The agreement is again very good and we 
obtain from the fit / = 0.46 (in the unit discussed in the introduction). As in the previous 
case we did the same comparison for 5 particles, obtaining again a very good agreement. 
Moreover also in this case the value of / is very close to that obtained for N = 2. Finally 
it is interesting to check if this agreement remains when N — > oo, i.e. for the stochastic 
irreversible equation (3.9). As can be seen from Fig. 7 the agreement is again very good 
and we still get the same value for the parameter / ~ 0.46. 

We were also able to compute ipk(r,E;2) and 4>k(r,E) for k = 2 and 3. It is also 
easy to compute the lowest order contribution to ipk(r,E;2) and 4>k(r, E), extending the 
computation from section 3. It is thus possible to compare, at least in this limited situation, 
the results. Contrary to what we found for k = and 1, ip2{f, E; 2) is quite different from 
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^{r, E;2). Analogously 4>2(r, E) and 02 ( r , E) differ significantly. A comparison of the 
term with k = 3 also shows deviations between the mechanical and the stochastic models 
although, surprisingly, much smaller than those found for k = 2. We note however that 
for this comparison we only have data for E = 0.012. 

4-3. The mean free flight. 

In kinetic theory one can define the mean free flight in two ways. Denoting by £i(X) the 
distance travelled by particle i before its first collision with an obstacle starting form the 
point X G iS/v, Zo is the average of £i(X) with respect to the SRB distribution (x + (dX, E; A) 
(it clearly does not depend on i). On the other hand we can consider the set S l N of points 
such that particle % is undergoing a collision, i.e. is on the boundary of one of the 
scatterers, then l\ is the average of £i(X) on S l N with respect to the projection of the SRB 
distribution /i + (dX, E; N). Observe that for the stochastic model these two quantities are 
identical. 

We computed both Z and h for the mechanical system with A = 2, 5, 50 and for the 
irreversible dynamics eq.(1.7) with E = 0.04. This was done by running a very long 
trajectory and taking the average of the distance travelled by a particle between two 
collisions to compute l\ or numerically integrating £i(X) along the trajectory to compute 
Zo- The results appears to be independent of A, at least within the accuracy of our 
computations, and are: 

Z = 0.46 
Zi = 0.58 

The value of Zo agrees very well with the value obtained from the fit of Z reported in 
the previous section. This implies that the correct way to compare the stochastic and 
the mechanical model is to use Zo as the mean free flight parameter in eq.(3.1). This is 
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consistent with the Green-Kubo formula eq.(1.3). We saw in sect. 2.1 that eq.(1.3) is 
well verified for the conductivity at small field of the deterministic model. In the case 
of the stochastic model eq.(1.3) reduces to an integral relation between F(°)(R, 0) and 
F<°)(R, 0*), see eq. (3.5) (3.6) in sect. 3. We did not prove this identity although numerical 
analysis for small iV seems to verify it. Finally the agreement between -i/> (r, 0; iV) and 
i/j (r,0;N) observed in sect. 4.1 tells us that the ratio between the conductivity for the 
deterministic and stochastic dynamics is independent of N at least for E — > 0. From 
eq.(3.7) we know that the conductivity for the stochastic model with one particle and 
E = is 3//4 so that also for the deterministic model we have 

«(0,1) = ^ (4-3) 

This relation is also very well verified by our computation for the 1 particle system. 

To better compare the deterministic and stochastic models we also computed the dis- 
tribution P(£, E; N) of £i(X) with respect to the SRB distribution. This distribution for 
5 particles and E = 0.04 is shown in Fig. 9 together with an exponential law with the 
same average, i.e. the distribution one would obtain running the same simulation for the 
stochastic case. We did similar computation for E = 0.04 and N = 2, 10 and 50. The 
results are again independent from N. 

5. Conclusions 

To put our study here in a physical context we note that a system of noninteracting 
electrons moving under the influence of an external electric field while undergoing elastic 
scatterings is often used as a crude model of electrical conduction in metals (the Drude 
model) [19], [20], [14]. To obtain the conductivity the velocity distribution function of the 
electrons is then computed from a Boltzmann type equation like eq.(3.2): with N = 1 and 
without the thermostatting E • J term. By doing this calculation only to linear order in 
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Fig. 9: Free path distribution P(l, 0.04; 5) compared with an exponential distribution 
with the same average 

E one avoids the problem that, without the thermostat eq.(3.2) does not have a solution 
since the system will never be in a true steady state [21]. A crucial ingredient in the 
calculation is the explicit assumption that for E = the distribution is one corresponding 
to equilibrium at a given specified temperature T, i.e. Maxwellian for a classical system. 

This description of the system of independent electrons interacting with the lattice of 
ions only via elastic collision is clearly not realistic. It is just used for obtaining a simple 
quick answer for the zero (small) field conductivity. For a more complete description of 
the steady state in a conductor one has to consider the system to be in contact with some 
reservoir which will absorb the heat generated by the current. It is this interaction with 
some external reservoir that was replaced, in the model considered here, by an artificial 
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thermostat. To our surprise however we found that this modeling does not lead to a 
Maxwellian distribution when E — > even when N is very large. This means that there is 
no equivalence of ensembles when it comes to modeling how the energy is extracted from 
the system- at least when there is no direct interactions between the particles other than 
that induced by the thermostat. We expect (and have some indication [22]) that this will 
change when we include collisions between the particles. Still it raises some caution about 
"thermostats" as a model for the description of stationary nonequilibrium states. 
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Figures' captions 

Fig. 1: General billiard structure with discs of radius R\ and Ri in a periodic box with 
side length 2L, N = 3 particles are shown. 

Fig. 2: Conductivity k(E, N) as a function of E for different N. 

Fig. 3: k(E,N) as a function of iV -1 for different E. Also plotted is the conductivity 
obtained from eq.(1.3) using the actual distribution function, see next section, for E = 0.04 
and compared with the value obtained by a direct simulation at the same field. Finally 
the highest line represents the conductivity obtained from eq.(1.3) using a micro canonical 
hypothesis. 

Fig. 4: Plot of 2Trripo(r : E; 2) for different values of E. The straight dashed line is 
obtained from the microcanonical distribution, Eq.(2.4). The dotted line gives the result 
for the stochastic model 

Fig. 5: Plot of 7r7"i/>i(r, E;2)/E for different values of E. The dotted line gives the result 
for the stochastic model 

Fig. 6: Plot of 27iripo(r, E; 50) for E = 0.04. Also shown are the results from simulations 
of (1.7) and from analytic solutions of the corresponding stochastic equation, Eq.(3.10). 
For comparison we also show the microcanonical result, corresponding to a Maxwellian. 
Fig. 7: Plot of 7rr^i(r, E; 50) /E and comparison with stochastic irreversible dynamics for 
E = 0.08 

Fig. 8: Comparison between the limiting value of the conductivity k(E) in the reversible 
model and in the irreversible model Koo(E). 

Fig. 9: Free path distribution P(l, 0.04; 5) compared with an exponential distribution 
with the same average 
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Fig. 6 
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Fig. 7 
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Fig. 8 
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Fig. 9 
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